function [Q,theta1Opt,output] = step1SR(theta1,setup,optim)
% This function implements the first (and thrid) step of the SR approach
% where all the Q parameters are found
setup.optimizer  = optim.optimizer;
theta1Values = struc2values(theta1,fieldnames(theta1));
if ~isfield(optim,'display')
    optim.display = 1;
end
if optim.optimizer == 1
    sigma              = optim.sigma;                 %The step size
    opts.SigmaMax      = optim.sigmaMax;              %The maximal value for sigma
    opts.LBounds       = setup.lowerBoundsValues;     %Lower bound for params
    opts.UBounds       = setup.upperBoundsValues;     %Upper bound for params
    opts.MaxIter       = optim.MaxIter*10;               %The maximum number of iterations
    opts.PopSize       = 3*length(theta1Values);      %The population size
    opts.VerboseModulo = 1;                           %Display results after every 10'th iteration
    opts.TolFun        = optim.TolFun/10;             %Function tolerance
    opts.TolX          = optim.TolX;                  %Tolerance in the parameters
    opts.Plotting      = 'off';                       %Dislpay plotting or not
    opts.Saving        = 'on';                        %Saving results
    opts.Resume        = 'off';
    opts.numCPUs       = setup.numCPUs;
    %theta1ValuesOpt    = cmaes_dsgeDisplay(@objectFunc,theta1Values,sigma,setup.InsigmaValues,opts,setup);
    theta1ValuesOpt    = cmaes_dsgeDisplay_mp_grendel(@objectFunc,theta1Values,sigma,setup.InsigmaValues,opts,setup);
elseif optim.optimizer == 2
    lambda0 = 10;
    theta1ValuesOpt = LMoptimizerTheta1(theta1Values,setup,lambda0,optim.TolFun/100,optim.TolX,optim.MaxIter,optim.display);
    %options = optimset('Display','iter','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX,'Jacobian','on');
    %theta1ValuesOpt = lsqnonlin(@nlsTheta1,theta1Values,setup.lowerBoundsValues,setup.upperBoundsValues,options,setup);
    
elseif optim.optimizer == 3
    if optim.display == 1
        options = optimset('Display','iter','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX);
    else
        options = optimset('Display','off','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX);
    end
    if isfield(optim,'numOptim')
        for i=1:optim.numOptim
            theta1ValuesOpt = fminsearch(@objectFunc,theta1Values,options,setup);
            theta1Values = theta1ValuesOpt;
        end
    else
        theta1ValuesOpt = fminsearch(@objectFunc,theta1Values,options,setup);
    end
elseif optim.optimizer == 4
    if isfield(optim,'numOptim')
        for i=1:optim.numOptim
            lambda0 = 10;
            theta1Valuestmp = LMoptimizerTheta1(theta1Values,setup,lambda0,optim.TolFun/10000,optim.TolX,optim.MaxIter,optim.display);
            %options = optimset('Display','iter','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX,'Jacobian','on');
            %theta1Valuestmp = lsqnonlin(@nlsTheta1,theta1Values,setup.lowerBoundsValues,setup.upperBoundsValues,options,setup);
            if optim.display == 1
                options = optimset('Display','iter','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX);
            else
                options = optimset('Display','off','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX);
            end
            theta1ValuesOpt = fminsearch(@objectFunc,theta1Valuestmp,options,setup);
            theta1Values = theta1ValuesOpt;
        end
    else
        lambda0 = 10;
        theta1Valuestmp = LMoptimizerTheta1(theta1Values,setup,lambda0,optim.TolFun/10000,optim.TolX,optim.MaxIter,optim.display);
        %options = optimset('Display','iter','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX,'Jacobian','on');
        %theta1Valuestmp = lsqnonlin(@nlsTheta1,theta1Values,setup.lowerBoundsValues,setup.upperBoundsValues,options,setup);
        if optim.display == 1
            options = optimset('Display','iter','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX);
        else
            options = optimset('Display','off','MaxIter',optim.MaxIter,'MaxFunEvals',optim.MaxFunEvals,'TolFun',optim.TolFun,'TolX',optim.TolX);
        end
        theta1ValuesOpt = fminsearch(@objectFunc,theta1Valuestmp,options,setup);
    end
else
    theta1ValuesOpt = theta1Values;
end
theta1Opt  = cell2struct(num2cell(theta1ValuesOpt),setup.selectTheta1(:));
[Q,output] = objectFunc(theta1ValuesOpt,setup);

end

